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Abstract 

Recent progress in the development of 
finite element methodology for the prediction of 
aerothermal loads is described. Two dimen- 
sional, inviscid computations are presented, but 
emphasis is placed on development of an approach 
extendable to three dimensional viscous flows. 
Research progress is described for: (1) utili- 

zation of a commercially available program to 
construct flow solution domains and display 
computational results, (2) development of an 
explicit Taylor-Galerkin solution algorithm, (3) 
closed form evaluation of finite element 
matrices, (4) vector computer programming 
strategies, and (5) validation of solutions. 

Two test problems of interest to NASA Langley 
aerothermal research are studied. Comparisons 
of finite element solutions for Mach 6 flow with 
other solution methods and experimental data 
validate fundamental capabilities of the 
approach for analyzing high speed inviscid 
compressible flows. 

Nomenclature 


A element area 

D solution domain 

Et total energy 

E,F flux components in x and y directions 

E,F artificial viscous flux components 

I J | determinant of Jacobian 

L length of element side 

M Mach number 

[M] mass matrix 

[N] element interpolation functions 

h unit normal vector to surface 

1 ,m components of unit normal vector 

p pressure 

r number of nodes per element 

S solution domain boundary 

s surface coordinate 

t time 

At time step 

u scalar variable, eq. (4) 

U typical conservation variable 

V element volume 

V velocity vector 

u,v velocity components in x and y 

di rections 

x,y,z coordinate directions 


Y ratio of specific heats 

p density 

£,n,c element natural coordinates 

Subscripts 

D constant element quantity, eq. (7) 

1 Index for node 

s constant surface quantity, eq. (10) 

« free stream 

Superscripts 

e element quantity 

n time step 


Introduction 


Thermal deformations and stresses Induced by 
aerodynamic heating on high speed vehicles are 
an important concern in structural design. 
Aerodynamic heating may have a significant 
effect on the performance of the structure, and 
effective techniques for predicting the heating 
and the thermal -structural responses are 
required. Typical flight vehicles that have 
significant aerothermal -structural interactions 
are shown in figure 1 and Include hypersonic 
cruise vehicles, the space shuttle, and advanced 
space transportation systems. Numerical methods 
such as finite difference methods and finite 
element methods play a significant role In 
fluid, thermal and structural analyses. Finite 
difference methods are the predominant numerical 
methods for computational fluid dynamics 
(CFD}*. A major problem in finite difference 
CFD is the generation of grids for realistic 
three dimensional vehicles such as finned pro- 
jectiles or complete aircraft. The development 
of well -constructed CFD grids has been the sub- 
ject of intensive research for several years. 
However, as noted in reference 2, the generation 
of a single grid that discretizes the entire 
flow region for a complex configuration is an 
extremely difficult and sometimes impossible 
task. These difficulties have led to recent 
research into alternative approaches for hand- 
ling complex geometries including zonal 

schemes^ and the finite volume approach^. 

Both finite differences and finite elements 
are used for thermal problems, but finite ele- 
ment methods are the predominant numerical 
methods for structural analysis^. Finite ele- 
ment methods are used extensively for structural 
analysis because of their capability to analyze 
a wide variety of realistic structures. The 
availability of general purpose finite element 
programs provide analysts a capability 
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to perform static and dynamic analyses with a 
spectrum of one, two, and three dimensional ele- 
ments that can be used to represent accurately 
virtually any structure. The integration of 
these analysis programs with modern computer 
graphics capabilities provides the basis for 
computer aided design. For these reasons, 
finite element methods are also finding increas- 
ingly wider acceptance for structural heat 
transfer analysis particularly where a thermal - 
stress analysis is of interest. 

Finite element methods are relatively 
immature for flow analysis in comparison to the 
finite difference method. Although finite ele- 
ments have been applied to flow problems since 
about 1974, major emphasis has been placed on 
incompressible flows. Reference 4 describes 
finite element approaches used for incompres- 
sible flows and provides the principal refer- 
ences to finite element flow literature through 
1982. A survey paper 6 summarizes progress of 
finite elements in fluid mechanics for the 
decade from 1974-1984. A 1982 assessment 6 of 
the finite element method's capability for the 
analysis of viscous compressible flows identi- 
fies early papers dealing in that area. Through 
1982 only modest efforts had been made to 
develop finite elements for high speed compres- 
sible flows for either inviscid or viscous 
fluids. Researchers in the Aerothermal Loads 
Branch at NASA Langley decided that finite ele- 
ment methods offered potential for further 
development particularly for applications to 
complex aerospace vehicles with strong aero- 
thermal -structural Interactions. A research 
program with this goal was initiated in 1983. 

The principal motivations for investigation of 
finite elements for compressible flows are its 
inherent capability to handle complex geometries 
and the basic ability to analyze fluid, thermal, 
and structural interactions with a single 
methodology. The researchers believe that 
finite element methods may help alleviate some 
of the difficulties encountered by finite 
difference methods in flow problems, but 
recognize that much research remains to be done 
before finite element methodology becomes 
competitive with finite difference methods for 
predicting aerodynamic heating on high speed 
flight vehicles. 

Research is underway at the NASA Langley 
Research Center to improve the capabilities and 
efficiency of finite element high speed compres- 
sible flow analysis methods and to develop more 
efficient integration of finite element fluid, 
thermal and structural analyses. The focus of 
the research is the prediction of aerothermal 
loads for complex three dimensional bodies. The 
research combines analyses with experimental 
studies conducted in the 8' High Temperature 
Tunnel (8' HTT) at NASA Langley 7 . The aero- 
thermal loads research includes aerodynamic 
heating on control surfaces^, pronounced loca- 
lized heating effects on wavy surfaces ^ and in 
shuttle tile gaps 1(1 . 

Finite element methodology for compressible 
flow is under development by a team of 
University, industry, and NASA researchers. 


The purpose of this paper is to describe some of 
the recent progress in the development of the 
finite element compressible flow analysis 
methodology. The paper details progress in two 
dimensional inviscid flow computations with 
emphasis on the approaches being investigated to 
provide efficient model generation and computa- 
tions for complex three dimensional flow 
domains. To develop the new approach, research 
is proceeding in five areas: (1) model genera- 

tion .and graphic displays, (2) algorithm devel- 
opment, (3) finite element matrices evaluations, 

(4) vector computer programming techniques, and 

(5) validation of the finite element approach. 
This paper highlights research in these areas. 

Model Generation and Display Capability 

The development of the model generation and 
display capability is based on currently avail- 
able methodology for finite element modeling of 
geometries. Finite element modeling of geomet- 
rical shapes is wel 1 -developed in structural 
mechanics, and this technology can be exploited 
to model and represent graphically the continuum 
models for complex flows. The commercially 
available PATRAN-G (Trademark of PDA Engineer- 
ing)* program 11 has been selected to generate 
flow models and display computational results. 
PATRAN offers mature state-of-the-art capability 
for geometrical modeling and results display for 
finite element analyses. Many of PATRAN's very 
general capabilities can beused advantageously 
for finite element flow analysis. 

Modeling solution domains with PATRAN begins 
with creation of a geometric representation of 
the computational domain. A two dimensional 
computational domain is represented by a number 
of regions, called patches, whose edges are 
described by parametric cubic lines. For a 
three dimensional computational domain, surface 
patches are used to establish solid regions 
called hyperpatches. A finite element mesh can 
then be created by dividing patches or hyper- 
patches into a number of nodes and elements. 
Different element types, such as triangles and 
quadrilaterals for 2D meshes and tetrahedrons 
and hexahedrons for 3D meshes can be assigned to 
suit the analysis and finite element program 
requirements. 

Examples of 2D and 3D computational domains 
are shown in figures 2 and 3. Figure 2a shows 
the patch representation of the flow field in a 
cove formed by a wing and elevon juncture. 

Figure 2b shows a very crude finite element mesh 
of the cove for illustrative purposes. Figure 
3a shows a cutaway view of a hyperpatch repre- 
sentation of the computational domain for flow 
around a cylindrical protuberance. The hyper- 
patch representation is shown with hidden lines 
removed and drawn with dotted lines, high- 


*Use of this commercial product does not imply 
endorsement by NASA. 
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lighting selected hyperpatches. Figure 3b shows 
a hidden-line view of the finite element mesh. 

Translator and inverse translator programs 
have been developed to interface PATRAN with the 
flow analysis. The translator program 12 uses 
the output file generated by PATRAN to produce 
the complete input data file required for a flow 
analysis. In addition to the basic input data 
of nodal coordinates and element nodal connec- 
tions, the translator generates the flow initial 
conditions, boundary conditions, and normal vec- 
tors along slip boundaries. This additional 
information is distinct from that required by 
finite element thermal -structural models cur- 
rently supported by PATRAN. After the flow 
analysis is completed, the solution is trans- 
lated into the format required by PATRAN for 
displaying computational results using the 
inverse translator program. Figure 4 shows 
schematically the computer usage for the flow 
analyses and the role of the translator and 
inverse translator programs. 


Finite Element Solution Algorithm 

Two general classes of explicit solution 
algorithms are being investigated; (1) Taylor- 
Galerkin algorithms, and (2) Petrov-Galerkin 
algorithms. The Taylor-Galerkin algorithm was 
first developed by Donea 13-14 and applied to 
compressible flows by Lohner, Morgan, and 
Zienkiewicz 15 "* 16 . Petrov-Galerkin algorithms 
have been applied to compressible Euler equa- 
tions by Hughes and Tezduyar 17 " 18 . A prelim- 
inary evaluation of these algorithms appears in 
reference 19. This paper utilizes the Taylor- 
Galerkin algorithm. 

The basic concept of Taylor-Galerkin 
algorithms is to use; (1) Taylor series expan- 
sions in time to establish recurrence relations 
for time marching, and (2) the method of 
weighted residuals with Galerkin's criteria to 
develop the finite element matrix equations 
describing the spatial distribution of the flow 
variables. A variant of a two time level 
Taylor-Galerkin algorithm proposed by Morgan and 
Lohner is used in this paper. 

The solution algorithm is applied to the 
basic flow equations in conservation form, 

3(11} . 3{E} , 3(F) _ n m 

S F~ + ST + ST ' 0 U) 


where p is density, u and v are the velocity 
components, and Et is the total energy. The 
pressure p is given by the ideal gas law, 

P - P (y- 1) CEt-1/2 (u2+v2)3 (3) 

Equation (1) is solved subject to appropriate 
initial and boundary conditions. Consider a 
solution domain D bounded by a surface S, (Fig. 
5.). The initial conditions consist of speci- 
fying the distributions of the conservation var- 
iables at time zero. Typical boundary condi- 
tions include: (1) specifying all conservation 

variables for supersonic Inflow on a segment of 
the boundary Sj, (2) requiring V»n ■ 0 on a 
slip boundary S2» where n Is a unit vector 
normal to the surface and V is the velocity vec- 
tor, and (3) using appropriate boundary condi- 
tions on the outflow surface S3. For super- 
sonic outflow surfaces, the conservation vari- 
ables are unknown and completely defined by the 
upstream flow. This natural boundary condition 
consists of surface integrals of flux components 
to be derived in the finite element formula- 
tion. For subsonic outflow surfaces there is a 
downstream dependency that can be satisfied by 
specifying pressure at the outflow plane. 


The solution domain 0 is divided into an 
arbitrary number of elements of r nodes each. 
Figure 5 shows typical quadrilateral elements 
(r=4) used in this paper. For simplicity, the 
finite element formulation will be given for a 
single scalar equation, 


3u , 3E 3£ 

3t 3x 3y 


= 0 


(4) 


where the variables u, E, and F are analogous to 
the corresponding vector quantities in eq. (1). 
Let (u} n denote the element nodal values of 
the flow variables u(x,y,t) at time t n . The 
time step At spans two typical times t n and 
t n+ i in the transient response. The computa- 
tion proceeds through two time levels, t«+i/2 
and t n+ i. At time level t n+ i/2* values for 
u that are constant within each element are com- 
puted explicitly. At time level t n+ j, the 
constant element values computed at the first 
time level are used to compute nodal values for 
u. In the time level t n +\ computations, ele- 
ment contributions are assembled to yield the 
global equations for nodal unknowns. The 
resulting equations are approximately diagona- 
lized to yield an explicit algorithm. 

Time Level t n +i/2 


where (U) is the vector of conservation var- 
iables, {E} and {F} are vectors of the flux com- To advance the solution to time level 

ponents of mass, momentum and energy. t n+I/2* a truncated Taylor series yields 


The vectors {U}, {E}, and (F), are given by 



p 


PU 


pv 


PU 


PU 2 +p 


pvu 

{U} = < 

pv 

tE}=< 

puv 

> , (F>=- 

2 

pv +p 


pE t 


(pE t +p)u 


(pE t +p)v 


» 4 


. j 


J 


( 2 ) 


u(x,y,t n+1/2 ) = u(x,y,t n ) 

+ j At|£(x,y,t n ) 


(5) 


Then eq. (4) is Introduced on the right hand 
side of eq. (5) so that 
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u(x,y f t n+ 1/2 ) = u(x,y,t n ) 

* % < S > 

At time level t n+ i/2t the dependent variable 
u(x,y,t n+ i /?) is assumed to have a constant 
value within an element. Since the 

dependent variable has the constant value 
u D n+1 ^ 2 within an element, the flux compo- 
nents E n+ */ 2 and F n+ 1/ 2 are constant within 
an element as well. 

At time level t n in the response u, E, and 
F vary within an element and are interpolated 
from nodal values. Thus, the following spatial 
approximations are used within an element. 


In advancing the solution to the next time 
level, the values of the dependent variables on 
the outflow surface will be required also. Let 
U s n+ l/2 denote the surface value on a 
typical element edge IJ on the surface S3 in 
figure 5. Following the approach used previ- 
ously, u s n+ l/ 2 is assumed constant on edge 
IJ at time t n+ i/2* but at time t n u, E and F 
vary along the edge. Thus the following approx 
imations are used on an element edge on S3, 


u e (x V t ) = u n+1/2 

u ' x, ” ,c n+l/2' u s 

(10a) 

u e (x,y,t n ) = [N(s )] {u} n 

(10b) 

E e (x,y,t n ) = CN(s)3 (E}" 

(10c) 

F e (x,y,t n ) = [N(s)3 {F} n 

(10d) 


u e (x,y,t n+ i/2) = u D n 1/2 

(7a) 

u e (x,y,t n ) * [N(x,y)] {u) n 

(7b) 

E e (x,y,t n ) = [N(x,y)] {E) n 

(7c) 

F e (x,y,t n ) = [N(x,y)] {F} n 

(7d) 


where [N(x,y)] denotes element interpolation 
functions and {u} n Is a vector of the element 
nodal quantities. For a typical quadrilateral 
element the interpolation functions are 
bilinear, and in a local natural coordinate 
system (see ref. 4) have the form 

N.= 1/4 (1+^5) ( l+n i n) 1=1, 2,3,4 (8) 

where and nj denote the nodal coordinates 
(5i, T )f*+l) of the element in the £-n plane. 

The equations for UQ n+ ^/ 2 for each element 
are derived by the method of weighted resi- 
duals. The spatial approximations given in 
eqs. (7) are introduced into eq. (6) to give a 
residual; the result is multiplied by a weight- 
ing function which in this case is unity. 
Finally, the weighted residual is integrated 
over the area A of the element. The result is, 

A u D n+1/2 = j [N] dA {u) n 

- r /, (SI « <»" ‘ 9 ’ 

- £ / A If “ > r '" 


With eq. (9), the dependent variable uq n+1 / 2 
for each element can be computed explicitly 
using nodal values for u, E and F from the prev- 
ious time t n . The constant element variable 
U g n +I/2 may be interpreted as a weighted 
average of an element's nodal values at time 
t n . A later section will discuss the evalua- 
tion of the three integrals that appear on the 
right hand side of eq. (9). 


where [N(s)] denotes the interpolation functions 
along an element edge. Using the method of 
weighted residuals, the values for u s n+1 /2 
are derived by integrating over the length L of 
an element edge. Hence 



Thus, eqs. (9) and (11) can be used to advance 
explicity the element and surface values of the 
dependent variables to t n +2/2* Beginning with 
nodal values of {u) n , {E}", and { F } n at 
time t n , eq. (9) is used to compute constant 
values u[) n+ l/2 for each element. In a simi- 
lar way, eq. (11) is used to compute constant 
surface values u s n+ l/ 2 for element edges on 
the outflow boundaries. These values are com- 
puted explicity by looping through all elements 
and element edges on outflow boundaries. 

Time Level t n +i 


To advance the solution to t n +j, forward 
and backward truncated Taylor series expansions 
at t n+ j/2 are used to write the approximation 

u (x.y,t n+1 ) = u(x,y,t n ) + At fjj. (x,y,t n+1/2 ) 

( 12 ) 


Then following the approach used previously, 
eq. (4) is introduced on the right hand side to 
yield 


u(x,y,t n+ i) = u(x,y,t n ) - At [ff(x,y,t n+1/2 ) 


+ |£( x »y»t n+ i/2) 


(13) 


At time levels t n and t n+ i the variables u, 

E and F are interpolated within an element by 
eqs. (7b - 7d), respectively. The equations for 
the nodal values for {u) n+1 can next be 
derived by the method of weighted residuals 
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using the interpolation functions N^(x,y) as 
weighting functions in the standard way.^ In 
this process, the last term in eq. (13) is 
integrated by parts. These operations yield 
equations for the nodal values of a single 
element. 


dA E <W> 


dA F(t n+1/2 ) (14) 


where 



CM] = 

/ [N] T [N] dA 

(15) 


'a 


and . 



(R} n+1/2 = 

-it / [l Ej +1/2 + m f" +1/2 ] 

(N) ds 
(16) 


[M] (U) 


n+1 


CM] {u} n 

“ /. I- 

( laN 

K fa 


+ At 


+ (R1 


n + 1/2 


In eq, (16) 1 and m are the components of the 
unit normal vector h. Following usual finite 
element procedures, the element matrices given 
in eqs, (14) and (16) can be assembled to form 
system (global) equations. 


surface, S 3 . The dependent variables 
u s n+ l/2 given by eq. (11) are used to com- 
pute Ee n *l and F s n+ 1 /2, The vector 
{R) n+ 1 / 2 represents the natural boundary con- 
dition for outflow surfaces mentioned previ- 
ously. No other boundary condition is imposed 
on supersonic outflow surfaces. 


The Taylor-Galerkin algorithm presented in 
the foregoing is based on the work of Lohner, 
Morgan, and Zienkiewicz^-lS. The principal 
differences in the present form of the algorithm 
are the introduction of the surface quantities 
represented by eq. ( 11 ) and the use of quadri- 
lateral elements. In reference 15 and 16, the 
boundary terms are treated differently and tri- 
angular elements are used. These references 
discuss similarities of the algorithm with 
second order finite difference methods such as 
Lax-Wendroff and a two step difference method by 
Burstein. The two time level Taylor-Galerkin 
algorithm is conditionally stable. For the 
diagonal mass matrix approach used herein the 
time step must satisfy the CFL criterion. 

For problems with strong shocks artificial 
viscosity is used to reduce oscillations typical 
of second order accurate approximations. The 
form used in this formulation is due to 
Lapidus^O, and consists of introducing arti- 
ficial viscous flux components E and F where 


E 


XA 


3U 

9X 


3U 

9X 


(18a) 


The matrix [M] defined by eq. (15) is cus- 
tomarily called the element consistent mass 
matrix. The term mass matrix is used because 
the matrix arises from the unsteady partial 
derivative 9u/9t which represents acceleration 
in structural dynamics. In structural dynamics, 
the coefficient matrix of the acceleration term 
actually has units of mass and is called the 
mass matrix. The term, consistent, is used to 
distinguish these matrices from diagonal mass 
matrices that arise from other discretization 
methods. A consistent mass matrix has off- 
diagonal terms that couple the element variables 
on the left hand side of eq. ( 14 )resul ting in an 
implicit scheme. In non-fluid applications, 
consistent mass matrices often result in 
improved time accuracy of solutions, although 
this observation has not been demonstrated 
conclusively for compressible fluids. 

To yield an explicit algorithm the consis- 
tent mass matrix [M] has been diagonalized for 
the present computations. The diagonal elements 
are computed from^ 

M. * / N, dA i = 1,2,... ,r (17) 

1 'A 1 

In the applications to be presented, the focus 
is on steady-state solutions so that time 
accuracy of the response computation Is not an 
i ssue. 

The last three terms on the right hand side 
of eq. (14) depend on element and surface values 
of the flux components evaluated at time level 
t n+l/2* Note that the last term {R} n+ 1 '2» 
given by eq. (16), uses values of the fluxes E 
and F on the element edges lying on the outflow 


F * XA 


9v 

ay 


9U 

ay 


U 8 b) 


where X is an adjustable constant, A is the ele- 
ment area, and U is a typical conservation vari- 
able. Computational experience has shown that 
1<X<2 provides acceptable solution convergence 
without excessive smearing of shocks. In the 
computations, the nodal values (U } n+1 at 
time level t n+ i are corrected (or smoothed) by 
incremental values { AU} computed from 


where 

[M]{AU> = 

■ (R) 



(19a) 

(R) = At / 
J \ 

' f en +1 (an. 

^ l j 9x 

+ F n+1 

| 9N 
jay 

] dA 

(19b) 


where E n+ 1 and F n+ 1 are the artificial vis- 
cous flux components evaluated at time level 
t n +i* The element integrals in eq. (19b) are 
evaluated using four point Gauss quadrature. 

Explicit Evaluation of Element Integrals 

Element integrals for the Taylor-Galerkin 
algorithm shown in eqs. (9), (11), and (14-16) 
were evaluated in closed form to avoid expensive 
numerical Integrations that are customary for 
quadrilateral elements. For two dimensional 
problems, these element integrals are either In 
the form of integration over the element area or 
along the element edge. As an example, the ele- 
ment mass matrix for a quadrilateral element, 
eq. (15), is given by 

[M] ■a [N] T [N] I J I d£ dn (20) 
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where [N] T is r the vector of the element inter- 
polation functions in terms of the element 
natural coordinates 5 and n. The determinant of 
the Jacobian jJ that appears in the above equa- 
tion represents the transformation from the ele- 
ment global coordinates x-y to natural coordi- 
nates £-n. The transformation permits the ele- 
ment integration to be evaluated over a 
square. The determinant of the Jacobian is 
given by 2 


3x 3y 
35 35 

3x dy_ 
3n 3n 


( 21 ) 


For a typical 4-node quadrilateral element shown 
in figure 5, the determinant of the Jacobian is 


k k 


j - 




,J,V. 


( 22 ) 


The difficulty in evaluating the mass mat- 
rix in closed form arises because the integrand 
involves the determinant of the Jacobian above, 
which if completely expanded produces an alge- 
braic expression with 128 terms. However, by 
realizing that because of the form of the inter- 
polation functions the determinant of the 
Jacobian varies with the natural coordinates £ 
and n in the form of second order polynomials, 
the mass matrix can be evaluated effectively by 
writing the determinant of the Jacobian in the 
form 


M=ZN.|J|. (23) 

Where |j|-j represents the value of the deter- 
minant of the Jacobian at the i^ node. 

Using eq. (8), the equations of |jh can be 
evaluated directly. For example, for a typical 
element node I (Fig. 5), 

|j|l = |J(t'-l,n-l)| 

= [(Xj-XjMyL-yjMxL-XjHyj-yj)]^ 

(24) 


Using the determinant of the Jacobian in the 
form of eq. (23), the mass matrix shown in eq. 
(20) can be evaluated in closed form. 

The approach described above can be extend- 
ed directly for the evaluation of three dimen- 
sional hexahedral element Integrals. Typical 
element integrals derived using the Taylor- 
Galerkin algorithm are in the same form as 
obtained for the quadri lateral element where 
the element integrations are either performed 
over the element volume or the element surface 
areas. As an example, the mass matrix for a 
hexahedral element is given by 



= f 1 J J £N] T END IjI dS dn dc (25) 
-1 - 1-1 11 

where 5,n,C are the element natural coordi- 
nates. The determinant of the Jacobian above, 
if completely expanded, contains approximately 
200,000 terms. To evaluate the integral in 
closed form, the determinant of the Jacobian 
was written in the form of eq. (23). Lagrange 
polynomials for a 27-node hexahedral element are 
used because they contain complete polynomials 
of 5*n,5 as they appear in the determinant of 
the Jacobian. 

The CPU time used by the closed form solu- 
tion of the mass matrix has been investigated 
and compared with CPU times required for differ- 
ent orders of Gauss numerical integration for 
quadrilateral and hexahedral elements. Although 
CPU time savings for quadrilateral elements are 
large (a factor of 50 for all integrals), the 
CPU time savings for hexahedral elements is even 
more significant. Typical CPU times for evalua- 
tion of a hexahedral element mass matrix are 
compared in figure 6. The figure shows that the 
closed form solution reduces CPU time for an 
element mass matrix significantly. Time savings 
in excess of an order of magnitude are obtained 
from the closed form solution in comparison with 
the popular 8 point Gauss integration method. 
Computational savings from explicit evaluation 
of other hexahedral element matrices are shown 
in figure 7. 

The time savings gained by explicit evalua- 
tion of the finite element integrals is an 
important step toward developing efficient 
finite element computations for large scale 
three-dimensional flow computations. 


Vector Programming Strategies 

The algorithms under evaluation are being 
implemented with vectori zation strategies speci- 
fically for the Langley VPS 32 (a Cyber 205 with 
16 million words of central memory). This com- 
puter achieves high computational speed when 
performing operations on long vectors. Vector 
lengths of at least 60 are required to justify 
vectorization efforts with maximum payoff 
achieved for vector lengths of 1000 or more. 

The predominant vector lengths in the vectoriza- 
tion scheme presented here are the number of 
elements in the finite element model with occa- 
sional operations using vector lengths equal to 
the number of nodes. 

The critical vectorization tasks are for 
those operations that are repetitive and 
performed at every time step. For finite ele- 
ment algorithms, these operations are: (1) 
assembly of element contributions from the right 
hand side of eq. (14) into the global system of 
equations, (2) solution of the global system of 
equations, and (3) application of boundary con- 
ditions. 
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The assembly of element contributions into 
the global system of equations is the process in 
finite elements which differs most from other 
numerical techniques and requires special 
routines for vectorization. Nodal unknowns are 
stored In one dimensional arrays from 1 to the 
number of nodes in the model, and in general, 
node numbering may be arbitrary throughout the 
mesh. Assembly of element contributions is per- 
formed using the VPS 32 FORTRAN-suppl led scatter 
routine which places an element contribution 
into the proper location in the system of equa- 
tions based on the element connectivity (node 
numbers of I,J,K,L in figure 5). Every element 
that contains a particular node in its connec- 
tivity provides Its own contribution to the sys- 
tem equations, therefore the assembly is an 
additive operation. "Scattering" alone would 
merely overwrite a previous element contribu- 
tion. The special vector routine, then does an 
"additive scatter." 

For an explicit scheme, solutions are 
obtained directly, so that operation (2) vector- 
izes naturally. Operation (3), the application 
of boundary conditions (fixed and slip) is an 
intrinsically scalar operation and difficult to 
vectorize. However, use of bit vectors to flag 
boundary nodes and use of VPS 32 FORTRAN 
supplied routines enables full vectorization of 
this operation. 

The program flowchart for the Taylor- 
Galerkin algorithm is shown in fig. 8. Using 
closed form integration and explicit vectoriza- 
tion, the program averages computational speeds 
of 4xl0“ 5 CPU sec/time step/node. The compu- 
tational savings derived from explicit vectori- 
zation of the boundary condition subroutine Is 
illustrated in Table I. 

Table I. - Computational Savings by 

Vectorization of Boundary Condition 
Routine 



UNVECTORIZED 

(SECS) 

VECTORIZED 

(SECS) 

CPU TIME PER 
CALL 

0.0066 

0.00058 

TOTAL CPU TIME 
1200 CALLS 

8.0 

0.7 

TOTAL PROGRAM 
CPU TIME 

36.0 

28.8 

PERCENT OF CPU 
TIME SPENT 

22.0% 

2.0% 

Validation 

of Finite Element Approach 


To validate the finite element approach 
described in the previous section, two test 
problems are presented. These test problems are 
chosen for their relevance to the prediction of 


aerodynamic heating and because of the 
availability of comparison solutions or 
experimental data. 

Panel Holder Leading Edge 

The first problem is Mach 6.57 flow over 
the blunt leading edge of the panel holder 22 , 
figure 9, used in the 8 1 HTT for testing a 
variety of configurations. The finite element 
model for the blunt leading edge of the panel 
holder at a 5° angle of attack, shown In figure 
10, consists of 3200 quadrilateral elements and 
3321 nodes. With the initial conditions of 
uniform Mach 6.57 flow and a maximum allowable 
time step of 0.5 x 10~ 7 seconds, a steady 
state solution was obtained in 3250 time steps. 
Convergence Is considered achieved when the I2 
norm of the change in density decreases 3 orders 
of magnitude. 

The density contours In figure 11 show the 
bow shock which forms ahead of the body and the 
expansion of the flow around the nose. Slight 
oscillations occur in the solution near the 
shock which is captured over 8 elements. 

For comparison, solutions were also 
obtained from two other numerical techniques. A 
method of lines 23 shock fitting solution was 
obtained for a cruder mesh containing one- 
seventh the number of nodes. An upwinded finite 
volume 24 shock capturing solution was obtained 
for a mesh similar to the finite element mesh. 
Comparison of these solutions with the finite 
element solution is shown in figures 12 and 13. 

A comparison of the velocity distributions along 
the mesh exit is shown In figure 12. The shock 
fitting technique adjusts the mesh boundary to 
coincide with the shock so that the solution is 
obtained only up to the shock, while the shock 
capturing solutions extend through the shock 
into the freestream region. All three solutions 
agree well. The oscillations in the finite ele- 
ment solution in the freestream region do not 
appear in the finite volume solution because it 
uses a more sophisticated smoothing technique. 
Also, no attempt was made to fine tune the 
finite element solution by increasing the arti- 
ficial viscosity. Figure 13 shows the distribu- 
tions of Mach number and pressure on the body. 
The finite element solution for Mach number 
falls between the shock fitting and finite 
volume solutions. Both the finite element and 
finite volume techniques predict a lower Mach 
number after the expansion point than the shock 
fitting solution. The pressure distributions 
are virtually identical for all three 
techniques. 

Wing-Elevon Cove 

The configuration for a wing-elevon cove is 
depicted In figure 14 which shows a sketch of 
the Space Shuttle Orbiter and a simplified 
cross-sectional view of its structure at the 
juncture between the wing and eleven. Without 
the spring-loaded rub seal, the pressure differ- 
ential across the windward and leeward surfaces 
of the wing would drive a portion of the hot 
boundary-layer flow into the cove where it would 
contact the aluminum load-bearing structure. To 
provide insight Into the problem and provide a 
data base for future vehicle design, a full- 
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scale model of the wing-el evon cove was tested 
in the 8' HTT. 25 Supporting analytical 
studies were also undertaken using the 2D 
compressible Navler-Stokes equations 8 , but 
only qualitative results were obtained. A 
partial cross section of the experimental model 
and instrumentation is shown in figure 15. The 
experimental parameters were the freestream 
conditions, the elevon deflection angle, and the 
cove leak area. For the finite element analysis 
a laminar flow case, test 11, in reference 25 
was chosen. The test parameters for this case 
are the freestream conditions shown In figure 
15, an elevon deflection angle of 25°, and a 
cove leak area equal to half the cove inlet 
area. The portion of the external hypersonic 
flow region used for the analysis is also shown 
in figure 15. 

The finite element model of the cove, shown 
in figure 16, consists of 1665 quadrilateral 
elements and 1802 nodes. The specified inflow 
conditions and the fixed cove exit pressure were 
calculated from experimental data from instru- 
mentation on the wing just ahead of the cove and 
just after the leak (figure 15). The Initial 
conditions in the hypersonic region were the 
inflow conditions and in the cove were cal- 
culated from experimental data at the cove 
inlet. The converged solution required 40,000 
time steps of 0.15E-6 seconds for a total CPU 
time of 3000 seconds on the VPS 32. A large 
number of time steps were required because the 
subsonic flow in the cove requires considerably 
more time to converge than the external super- 
sonic flow. 

The finite element solution is shown in 
figures 17 and 18. The Mach number contours 
(figure 17) show oscillations in the solution 
ahead of the shock and in the expansion region 
at the cove inlet. The solution shows the 
expansion of the flow around the wing trailing 
edge. A portion of the external flow enters the 
cove and flows subsonically (M^0.4) through the 
leaking seal. The bulk of the external flow 
impinges on the elevon forming an oblique 
shock. The pressure contours (figure 18) are 
relatively smooth and clearly define the oblique 
shock which was captured over 8 elements. 

Figure 19 compares the finite element solution 
for pressure along the inner cove wall and 
elevon with experimental data. The experimental 
data shows that the pressure is nearly constant 
in the cove and rises gradually along the 
elevon. The finite element solution predicts 
the constant cove pressure with a steep pressure 
rise along the elevon. Figure 19 also shows the 
pressure on the elevon as approximated by 
oblique shock theory for a wedge angle of 25°. 
Downstream of the computational domain, the 
experimental data appears to be asymptotical ly 
approaching the oblique shock theory pressure. 
The finite element solution slightly over- 
estimates the oblique shock theory pressure. 

The actual flow phenomena is complicated at the 
cove entrance by viscous effects not included in 
the present analysis. Viscous effects may 
account for the disagreement In the pressure 
predicted by the inviscid finite element 
analysis and the experimental data along the 
elevon (s>7). 


Concluding Remarks 

Recent progress in the development of 
finite element compressible flow methodology is 
described. The goal of the research is to 
provide an effective computational procedure for 
complex flight vehicles with significant 
aerothermal -structural interactions. Two 
dimensional inviscid flow computations are 
emphasized, and the approaches being investi- 
gated to provide efficient mesh generation and 
effective computations extendable to three 
dimensional flows are described. The mesh gen- 
eration approach is based on currently existing 
finite element methodology used in structural 
mechanics. The commerically available PATRAN-G 
program is used, and many of PATRAN’s very gen- 
eral capabilities for constructing representa- 
tions of solution domains can be used advantage- 
ously for finite element flow analysis. To 
interface PATRAN with the flow analysis pro- 
grams, interface programs have been developed. 

A two-step, explicit Tayl or-Galerkin solu- 
tion algorithm is described. The algorithm is 
similar to earlier second order finite differ- 
ence schemes and uses artificial viscosity to 
smooth oscillations in the vicinity of shocks. 
The element integrals that appear in the 
Taylor-Galerkin algorithm are evaluated in 
closed form for two dimensional quadrilateral 
and three dimensional hexahedral elements. 
Numerical calcuations show that CPU times can be 
reduced significantly by the closed form 
evaluation of the matrices. The algorithm has 
been implemented on the Langley VPS 32 vector 
computer with special programming strategies to 
insure computations with long vectors. Typical 
vector lengths equal the number of nodes or 
elements are employed. 

To validate the approach, two test problems 
of interest to aerothermal research at NASA 
Langley were studied. The first problem is Mach 
6.57 flow over the leading edge of a panel 
holder used in the NASA Langley 8' High Tempera- 
ture Tunnel (8' HTT). The finite element solu- 
tion is compared with a method of lines shock 
fitting solution and an upwinded finite volume 
shock capturing solution. The finite element 
solution smears the shock over about 8 elements 
and shows slight oscillations. The solutions 
generally show good agreement for Mach number 
and pressure on the body surface and velocity 
components on the outflow plane. 

The second test problem considers a leaking 
space shuttle wing elevon cove. The effects of 
a leaking wing elevon seal have been studied 
previously in the 8* HTT using a full scale 
model. The complex cove geometry combined with 
Mach 6 external flow and subsonic cove flow have 
taxed previous computational approaches to the 
problem. The finite element solution modeled 
the complete, highly irregular cove geometry and 
captured essential features of the flow. The 
external shock on the wing was captured and the 
low Mach number subsonic flow in the cove was 
represented qualitatively. Comparison of the 
finite element predicted pressure distribution 
along the cove inner wall shows excellent agree- 
ment with experimental data within the cove, but 
the finite element solution overestimates the 
pressures on the elevon surface in the 
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supersonic outer flow. A possible source of 
this discrepancy is viscous flow effects at the 
cove entrance that are not included in the 
inviscid finite element analysis. 

The two test problems have validated 
fundamental capabilities of the finite element 
approach for analyzing high speed inviscid 
compressible flows. The ground work has been 
established for further investigation of the 
approach for three dimensional flows and the 
extension to include viscous effects. Much 
research remains to be done to establish a 
robust, generally effective method for complex 
flight vehicles with strong aerothermal -struc- 
tural interactions, but progress towards achiev- 
ing this goal has been demonstrated. 
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Fig. 1 Flight vehicles with aerothermal -structural interactions. 



a) Patch representation of computa- 
tional domain 




a) Hyperpatch representation of 
computational domain 



COMPUTER USAGE FOR FLOW ANALYSIS 



Fig. 4 Computer usage for flow analysis. 
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Fig. 5 Two dimensional finite element 
computational domain. 



Fig. 6 Comparison of execution times for 
evaluation of hexahedral element 
mass matrix. 


8 Gauss pts 





Fig. 8 Program flowchart for Taylor- 
Galerkin algorithm. 



Fig. 9 Panel holder in 8- Foot High 
Temperature Tunnel (8 * HTT). 


Fig. 7 Comparison of execution times for 
evaluation of hexahedral element 
matrices. 
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3321 nodes 

3200 quadrilateral elements 



4R 


r — R — H 

Fig. 10 Finite element mesh for blunt 

leading edge of 8' HTT panel holder. 



Fig. 12 Comparative solutions at outflow plane 
for flow over blunt leading edge of 
8' HTT panel holder. 



Fig. 11 Finite element predicted density 

contours for flow over blunt leading 
edge of 8‘ HTT panel holder 

(units are xlO" 6 lb m /1n 3 ). 



Fig. 13 Comparative solutions along the body surface 
for flow over blunt leading edge of 8* HTT 
panel holder. 
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Fig. 14 Space Shuttle wing-elevon cove 



Fig. 15 Partial cross-section of experimental cove 
model . 



Fig. 16 Finite element mesh for flow in wing 
elevon cove. 
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